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Multiple time correlation functions are found in the dynamical description of different phenomena. 
They encode and describe the fluctuations of the dynamical variables of a system. In this paper 
we formulate a theory of non-Markovian multiple-time correlation functions (MTCF) for a wide 
class of systems. We derive the dynamical equation of the reduced propagator, an object that evolve 
state vectors of the system conditioned to the dynamics of its environment, which is not necessarily 
■^j- ' at the vacuum state at the initial time. Such reduced propagator is the essential piece to obtain 

i multiple-time correlation functions. An average over the different environmental histories of the 

f^*) . reduced propagator permits us to obtain the evolution equations of the multiple-time correlation 

£SJ ' functions. We also study the evolution of MTCF within the weak coupling limit and it is shown that 

_ | the multiple-time correlation function of some observables satisfy the Quantum Regression Theorem 

(QRT), whereas other correlations do not. We set the conditions under which the correlations satisfy 
^fy • the QRT. We illustrate the theory in two different cases; first, solving an exact model for which 

the MTCF are explicitly given, and second, presenting the results of a numerical integration for a 
, system coupled with a dissipative environment through a non-diagonal interaction. 
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Introduction and motivation. Many research contexts are focused on the dynamics of a system (S) that 
is affected by an environment (B) from which it cannot be considered isolated. Examples of such situations are 
encountered in statistical physics, condensed matter and quantum optics. We found a concrete example in the 
description of the dynamics of an atom (<S) immersed in an electromagnetic field (B) • 

In some circumstances, the analysis of the dynamics of the system is done using the expectation values of its 
observables over state vectors of the whole system, and then averaging over the environmental degrees of freedom. 
However, in some other situations, like when studying the response of a system to an external EM field, some additional 
fi || information is needed. In particular, for the analysis of the spectroscopic properties of a system some multiple-time 
+Jj ■ correlation function (MTCF) has to be computed, usually a two-time correlation function. 

The dynamics of the system S is usually described through its reduced density operator. Such operator verifies 
some master equation that in the Markovian case is of Lindblad type @, EL IE1S 0- Complementary to the 
q-( master equation approach, a series of Monte-Carlo type of approaches based on the so called stochastic Schrddinger 
equations 0> @, 0> ED have been developed in the last decade. In such schemes, the dynamics of system state 
, £^ ' vectors is integrated, and after an average is made over many realizations of environment histories that eventually 
are understood as a noise and takes into account the environment influence on S. In the non-Markovian case, within 
the context of nuclear magnetic resonance, the Redfield master equation was developed 0,0- On the other hand, 
several non-Markovian stochastic Scrodinger equations have been established and studied very recently for systems 
influenced by a structured environment 0, Q, 0, 0, 0, 0, 0|. In the Markovian case, an important tool to 
compute time correlation functions is the Quantum Regression Theorem (QRT) 0,0 El- The theory of stochastic 
Schrodinger equations , initially elaborated to compute the expectation values of system observables, has been extended 
by many groups [ij, l2ll \2'A |23| to calculate multiple-time correlation functions for the Markovian case. In addition, 
such stochastic methods agree with the results expected from the QRT. It is then natural to develop an equivalent 
theory within the stochastic Schrodinger equations approach, of multiple-time correlation functions for non-Markovian 
interactions. Interest arises particularly in those cases where the QRT is not valid. 

Along with the theoretical interest, there is a practical need for a theory of non-Markovian MTCFs, in order to 
apply them to systems with experimental interest. Because of its potential applications, Photonic Band Gap (PBG) 
materials constitute an interesting example, in which the electromagnetic field displays a band-gap structure |24ll25| . 
As a consequence, the correlation function of the field is highly non-Markovian, and therefore an atom in interaction 
with such field will display non-Markovian dynamics. In this context, there are several works devoted to investigate 
the dynamics of few-level atoms in PGB-materials Efil For them, non-Markovian stochastic Schrodinger 

equations have been successfully applied _28j. Therefore, the theory of MTCF we intend to derive in this letter is 
important from both experimental and theoretical points of view. 

Multiple-time correlation functions. A frequently used Hamiltonian model in the study of the dynamics of 
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S with Hamiltonian Hs, in interaction with B is 

H = H s + LB^ + L^B + H B = 

H S + X(L^2 9ndl + i f Sn a n J + ^ W ™ a « a ™> ( X ) 
^ n n ' n 

where the operator L acts on the Hilbert space of the system, a„,a^ arc the annihilation and creation operators on 
the environment Hilbert space and A is a suitable perturbation parameter that eventually can be taken equal to one. 
The g' n s are the coupling constants and the ui' n s are the frequencies of the harmonic oscillators that constitute the 
environment |2flj . 

We are interested in the evaluation of N-time correlation functions, defined for a set of observables 
{Ai(t±), • • • , A^^n)} = A(t) in Heisenberg representation as 

C A (t|*o) = (*o|A(ii)---^iv(*iv)|*o), (2) 

with t\ > t2 > ■ ■ ■ > tN and t = {t\, • • • , t]y}- The initial state of the full system is taken as the tensor product of a 
system state l^o) and the environment state \zq), i.e. \^o) — \ipo)\zo). 

In the partial interaction picture with respect to the environment, the N-time correlation function is defined as 
CA_(t\if?o) = (^ol IIt=i ^7~ 1 (^*' 0)^4jZ^j(tj , 0)|5'o)? where Uj is the evolution operator of the system in the interaction 
picture. A suitable basis to treat the environment B is a coherent state basis, \zi,z 2 ,--- ,z„,---) = \z) in the 
Bargmann representation P, l30| . In such basis the resolution of the identity is given by 1 = j d(j,(z)\z) (z\ with 
d/j,(z) = Y[iLi(d 2 Zi exp(— |zi| 2 )/7r), and when inserted in the definition of CaW^o) it follows 

N 

C A (t|* )= d^^olG-'^^YlAG^i + l)]^), (3) 

i=i 

with to = 0, tpj + i — and zn+i — zq. We have introduced the reduced propagators G(i,i + 1) = G(z? Zj+i|titj+i) = 
(zi\Ui(ti, ti+i)\zi+i), which act on the system Hilbert space, giving the evolution of system state vectors from U + \ 
to t{, conditioned that in the same time interval the environment coordinates go from Zi+i to Zj. Once their time 
evolution is solved, the time-correlation function can be obtained. Therefore, to proceed further we need to derive 
the equation of motion of the reduced propagator G(i,i + 1), by considering its time derivative with respect to ij. 
Taking into account the fact that the evolution operator Uj satisfies the Schrodinger equation in the partial interaction 
picture, and after some manipulations, we arrive to the equation 

dG{l ^ l) =(-*H s + Lz* t . LU i+hti )G(i, i + 1) 

-L 1 / dTOt(ti-T)(zi\Ui(ti,ti +1 )L(T,ti +1 )\z i+ t), (4) 



with L(t',t) = e iH B t e -iH{t-t') Le iH{t-t') e -iH B t' At = i Y Jn g nZi ^ e ^ an d a (t - t) = E„ l3«| 2 e" l "" ( ^ r) . The 
function z^t is a sum over time dependent coherent states and a(t — s) is its time autocorrelation function, as it can 
be verified by computing the average M[zijz* s ] with respect to the measure dfJ>(z). When integrating equation 
we still have the problem that it is not possible to compute exactly the matrix element (zi\Ui(ti, ^+i)L(r, ti+i)\zi+\) , 
and to express it as a function of the reduced propagator. This would bring to an explicit equation for the reduced 
propagator. Since only in very exceptional cases exact solutions can be obtained, some approximate scheme has to 
be taken. One possible way is to treat L(r, in the weak coupling limit. Else, sometimes it is possible to assume 
that (zi\Ui{ti, ti + i)L(r, ti + i)\zi + i) — 0{zi+\Zi, T)G(i, i + 1), where the operator O has to be constructed [3l]]. In 
this case we have 



gg(M+i) i .„ , T * rt 



—if J dra(U - T)0(z i+1 Zi,t i+1 ,T)^G(i,i + 1). (5) 

The equation (0J or its approximate versions, in particular equation (JjjJ, depends on two time dependent functions, 
z* t . and z* +lt ., that take into account the "history" of the environment and lead to a conditioned dynamics of 
the system with respect to the environment dynamics. They constitute one of the results of this letter and they 
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are the starting point to compute the MTCF in the non-Markovian case. The integration of the equations for the 
reduced propagators along with their initial conditions, G(i,i + 1) = exp(z*Zi+i) , leads to the evaluation of the 
iV-time correlation functions previously defined. In addition, let us comment that since the equation for the reduced 
propagator Q is made for an initial state of the environment different form the vacuum, this equation permits us to 
evaluate expectation values and correlation functions of system observables with more general initial conditions that 
the one usually taken, i.e. \^&o) — I "00 ) I vacuum), but we will not report those results here. 

In principle, and under the conditions we have considered, the non-Markovian multiple-time correlation functions 
derived in this section are rather general. We shall look closer to the time correlation functions we have derived, 
showing how they can be related to the Quantum Regression Theorem which applies for the Markovian case, and we 
will show a couple of examples. 

Beyond the Quantum Regression Theorem. Weak-coupling limit. Once we have the multiple time 
correlation functions, we may compute them directly from the stochastic method. Nonetheless, this may turn to be 
an expensive strategy from the numerical point of view, what is specially true when the number of environmental 
degrees of freedoms needed to correctly describe its correlation function is large. For those cases, we present in this 
section a set of coupled differential equations in which the stochastic average has been done analitically, and which 
evolve the non-Markovian two-time correlations up to second order in a convenient perturbative parameter A. 

The method we will follow consists in deriving the reduced two-time correlation (ipo 
G 1 ^ (ziO\t'0)V t > AG 1 (zlz 2 \t't)V t BGt,o{z%0\tO) | ip ) with respect to t', and then performing analitically the av- 
erage over the variables z\ and z%. In order to do that, it is necessary to use a perturbative expansion of the 
propagators in the interaction picture with respect to the system, which is used because it leads to more handable 
expressions. We then arrive to the following set of differential equations for the two-time correlation functions up to 
0(X 3 ) 

A(*o I A{t')B{t) | *„) = (V \(i{[H s ,A}}(t')B(t) 

+ f dTa*(t'-T){V T - t >tf{A,L}}(t')B(t) 
Jo 

+ f dra(t' -T){[tf,A]Vr- t ,L}(t')B(t) 
Jo 

+ fdra(t' -T){[L\A]}(t'){[B,V T -tL]}(t))\*o). (6) 

In this expression, the time dependencies of operators are now in the total interaction image evolution operator. We 
also denote V t >L = exp{i£si'}L = exp (iHst')Lexp (—iHsf), where V? = exp{iCst'} is the free system Liouville 
operator, acting in the two sides of the immediately contiguous system operator. It can be checked that only when the 
last term of 10 vanishes, i.e. when [X^A] = or [B,V T -tL] — 0, the Quantum Regression Theorem applies. Notice 
also that this term is zero in the Markovian case, since the corresponding correlation function a(t' — r) = T5(t' — r) 
vanishes in the domain of integration from to t. In summary, the previous equations leads to the computation of 
the MTCF and they contain the conditions under which the QRT remains valid in the weak coupling limit. Equation 
© is another result of this letter. 

A solvable example. To illustrate the theory proposed in this letter, we shall apply it to a simple solvable model 
described by the hamiltonian with L — a z and H$ — %o~z- This model describes the dynamics of system state 
vectors towards one of the eigenstates of the system Hamiltonian. Notice that in this case [H$, L] = and then O = L 
(see equation EJ . 

Let us consider the two-time correlation of A = {{0, a}, {(3, 0}} and B = {{1, 0}, {— 1, 0}} = a z . For an initial 
system state ji/'o) = vpoi I + ) + ^02 1 — ) the result is 

C AB (t't\* ) = e-^'^ip^oie-^' - aft^ 2 e" f '}, (?) 

with the definition Jj|!(a) = f^dr J d dsa(r — s). For the case in which A = B = <j z , we have C 0z a z = 1. Another 
type of two-time correlation function is the one corresponding to the observables A — {{0, a}, {(3, 0}} and B = 
{{0,a'},{/3',0}}givenby 

C AB {t't) = e^'^a/^oilV^ 4 '-') + a'Pl^fe-^'-^}, (8) 

with D(t't) = I^ia*) + l( t T (a) + Ifo(a) + I$f (a) - 7 t * *(a) - I$(ot). In the first case, since [L, B] = the last term 
of equation (JSJ is zero, and the QRT is valid; meanwhile, in the second type of correlation functions neither [L, B] = 
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nor [X^, A] = and therefore the QRT is not expected to hold. This can be verified by direct substitution of the 
two-time correlation functions obtained into JJJ, giving further support to the equations proposed. 
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FIG. 1: Figure (a) represents the imaginary part of C CTxC r 2 , the solid line represents the analytical result JJj, that in this case 
coincides the Quantum Regression Theorem result. Dotted an dashed lines represent the stochastic result with an average over 
10 2 and 10 5 trajectories respectively. Figure (b) represents the imaginary part of C CTajCT . Comparing the result of the QRT 
(dotted line), to the exact result given by © (solid line), it is clear that QRT is not valid in this case. The dot-dashed and the 
long-dashed lines are obtained averaging over 10 2 and 10 4 trajectories respectively. 

As an illustration, the hgure (JIJ shows two-time correlation functions of the system C CTx(Ti and C ax<r , with two 
oscillators in the environment with parameters 31,32 = .9 = 1 arid u>\ = 6, LO2 = 2. The initial system state taken in 
all computations in this letter is l^o) = |V'o)l v acuum) with |^o) = ((1 + 2i)|+) + (1 + )) /V7. It is clear from the 
figures (c) and (d) that the QRT does not apply for C as . ayl since [L, B] ^ and [L\ A\ ^ 0. 

An example of dissipative system. Let us now apply the theory derived in this paper for the case L = a i2 , 
and for a dissipative interaction. Within the perturbative approximation, the operator 0(z,t,r) can be replaced by 
its zero order perturbative expansion, V T -tL = u\2 exp {iu{t — t)}, where u> is the system rotating frequency. We 
propose the following correlation function, a(t — r) = $^m=— W2 C( m ) e ~ ITrm ( t_T )/ T ; with the coefficients, C(m) = 

2T I-t dta(t)e l7Tmt / T , which represents the Fourier series of the function (r/2) exp {T\t — r|}. In those equations, T 
is the time window in which the correlation function is expanded in the series. The more members we add in the 
former sum, the closer the solution is to the exponential decaying correlation function, and the larger we can fix the 
recurrence time T. The decaying behaviour is displayed in Figure (J2J), representingthc correlation C<j x(Tx as a function 
of t' — t = r. In this figure, we compare the evolution given by © to the ensemble averaged stochastic evolution for 
different number of trajectories. We also consider the correlation function a(t — r) = (r/2) exp {— T\t — r|}, in order 
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FIG. 2: Two time correlation C CTx(Tx (t' , t) for the coupling L = o\2, and the Fourier series of the exponential correlation function 
with v — 8 oscillators. The parameters are: u = 0.1, V = 1, perturbative parameter g = 0.2, recurrence time T = 40, and 
initial time for the correlation t = 1. Solid line represents the solution of the system ©, whereas dotted, dashed and long 
dashed lines gives respectively the result of the stochastic method for k = 5x 10 s , 10 6 , and 6 x 10 7 trajectories. An increasing 
accordance with the system curve (equation^ is observed as the number of trajectories grows. 

to study the validity of the QRT for Co- x(Tx (see Figure |3J). 

In summary, we derive in this letter a theory for non-Markovian multiple-time correlation functions from stochastic 
Schrodinger equations. The starting point of the theory is the evolution equation for the reduced propagator, that 
evolves vectors in the Hilbert space of the system S conditioned to the dynamics of the environment. Remarkably, 
such equation depends on two different histories of the bath. With the reduced propagator, the multiple time 
correlation functions are formally obtained. Furthermore we have derived, in the weak-coupling limit, the set of 
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FIG. 3: Two time correlation Cxx = Ca m a m (see text) The parameters are: ui = 0.1, r = 1, perturbative parameter g — 0.4, 
and initial time for the correlation t = 10. Solid line represents the solution of the system ©, and dotted line gives the result 
expected with the QRT. Because the last term in (|SJ is non zero, both results are different from each other, and the QRT is 
not valid. 

coupled differential equations that satisfy the two-time correlation functions. This set of equations is a generalisation 
of the Quantum Regression Theorem, and they show clearly the cases in which such theorem is not longer valid. We 
have verified the theory by applying it to two cases. For a solvable model, we computed the two-time correlation 
functions explicitly, displaying a case in which the QRT is fulfilled and a case in which it is not. For a case in 
which the system has a non-diagonal interaction with the bath, and the bath has a decaying exponential function, we 
have numerically integrated the multiple-time correlation functions. Following the procedures here established, higher 
order time correlation functions might be calculated. We believe that this work is relevant, aside from its intrinsic 
theoretical interest, in the description of the dynamics of small systems, such as atoms immersed in photonic crystals 
as well as other situations where non-Mar kovian effects are relevant. 
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